import numpy as np
import environment_pool as env

def benchmark(x,npar):
    return env.benchmark_pool(x,npar)
## approximately calculate the gradient through the difference method
# notice that if the benchmark is replaced with your own cost function, the
# input parameters of this function should be mapped.
def derivation(x,deltax,npar):
    y_n= np.zeros((npar))
    x_n=x.copy()
    #print('x_n:',x_n)
    for i in np.arange(npar):   
        x_n[i] = x[i] + deltax
        y_n[i],match_request,_ = benchmark(x_n,npar)
    y,match_request,_ = benchmark(x,npar)
    return (y_n - y) / (deltax)
def SMALL0(value):
    value2 = value.copy()
    value2[np.abs(value)<1e-15] = 1e-15
    return value2
# value = np.array([1e-20,0.3,0.2,-1e-21,0.4])
# print(value)
# print(np.log(SMALL0(value)+0.0j))
## the calculation of the current i3 in condition of time-varying field
#more information please refer to the referrence.
def solve_i(L2,dx,h,dfx,i2o,i2oo,i3o,i3oo,R1,R2,R3,B,S3):
    deltaS = h * dx + 0.5 * dfx * dx ** 2 # deltaS is approximately calculated by the area of a trapezoid.
    L3 = S3
    Z2 = np.abs(R2 + L2 * (np.log(SMALL0(i2o)+0.0j) - np.log(SMALL0(i2oo)+0.0j)))
    # the value of np.abs(1/i2*(di2/dt)) is transformed to be d(np.abs(np.log(i2)))/dt, 
    # and the gradient is approximately calculated through the difference.
    Z3 = np.abs(R3 + L3 * (np.log(SMALL0(i3o)+0.0j) - np.log(SMALL0(i3oo)+0.0j)))
    # print(Z2.shape,Z3.shape)
    Ztotal = R1 + Z2 * Z3 / (Z2 + Z3)
    i3n = (dfx / Ztotal * Z2 + B * deltaS) / (Z2 + Z3)
    i2n = dfx / Ztotal - i3n
    i3n = (i3n) / (np.abs(i3n) + np.abs(i2n)) #normalization
    i2n = (i2n) / (np.abs(i3n) + np.abs(i2n))
    # i3n = (i3n) / (max(np.abs(i3n),np.abs(i2n))) #normalization
    # i2n = (i2n) / (max(np.abs(i3n),np.abs(i2n)))
    # 
    return (i2n,i3n)
def Geodesic_cfg_and_run(root,data):
    '''
    data = {
    '<freq_MHz>':1268,
    '<file_name>':'results',
    '<mesh_file>':'A320-mesh.nas',  #"A320-mesh_new.facet"
    '<Rx>':-7.90782e3,
    '<Ry>':-0.585195e3,
    '<Rz>':-0.022963e3,
    '<Tx>':-11.642e3,
    '<Ty>':1.55778e3,
    '<Tz>':2.885e3,
}
results = Geodesic_cfg_and_run("C://Users//Administrator//Desktop//Geodesic",data)
print(results)
'''
    import os
    from template_cfg import GeodesicConfig
    xdsgen = GeodesicConfig('Geodesic_template')
    xdsgen.set_geodesic_root(root)
    xdsgen.set_config(data)
    # print(f'nowPath1:{os.getcwd()}')
    xdsgen.load_template('MEDO_template.json')
    nowPath = os.getcwd()
    temp_path = f"{root}//temp"
    import shutil
    # print(data["<mesh_file>"])
    shutil.copy(data["<mesh_file>"],f'{temp_path}//{data["<mesh_file>"]}')
    os.chdir(temp_path)
    xdsgen.set_cfg()
    xdsgen.write_cfg('C://temp.json')
    xdsgen.write_cfg('temp.json')
    xdsgen.run_meshgen()
    results = xdsgen.results()
    os.chdir(nowPath)
    # print(f'nowPath2:{os.getcwd()}')
    return results
def file2lines(file_path,encoding=None):
    f = open(file_path,encoding=encoding,errors='ignore')
    lines = f.readlines()
    f.close()
    return lines
def double_size(array,cur_size,grow_rate=0.2):
	if(np.size(array)==0):
		return array
	nsize = np.size(array[:,0])
	msize = np.size(array[0,:])
	if(cur_size >= nsize*0.95):
		nsize = int(max(cur_size,nsize*(1.0+grow_rate)))
		new_array = np.resize(array,(nsize,msize))
		new_array[cur_size:,]=0.0
		return new_array
	return array
mapNAS2NP = {
    'CBAR  ':[2,'segs'],
    'CTRIA3':[3,'tris'],
    'CTETRA':[4,'tets'],
    'CQUAD4':[4,'quads'],
    'CPENTA':[6,'prisms'],
    'CTRIA6':[6,'tri6s'],
}
def read_nas(nasPath):
    nasLines = file2lines((nasPath))
    mapResults = {}
    mapResults['pts_num'] = 10000
    for type in mapNAS2NP.keys():
        elems_strs = mapNAS2NP[type][1]
        vertex_num = mapNAS2NP[type][0]
        mapResults['n%s'%elems_strs] = 10000
        elems_num = mapResults['n%s'%elems_strs]
        mapResults[elems_strs] = np.zeros((elems_num,vertex_num+1),dtype=np.int64)
        mapResults['i%s'%elems_strs]=0
    mapResults['pts'] = np.zeros((mapResults['pts_num'],3))
    mapResults['ipt'] = 0
    max_ipt = 0
    for nasline in nasLines:
        if(nasline[0]=='$'):
            continue
        if(0==nasline.find('GRID*')):
            mapResults['ipt'] = int(nasline[5:24])
            iPt = mapResults['ipt']
            max_ipt = np.max([max_ipt,iPt])
            mapResults['pts'] = double_size(mapResults['pts'],iPt,1.0)
            mapResults['pts'][iPt-1,0]=float(nasline[26:56])
            mapResults['pts'][iPt-1,1]=float(nasline[56:73])
        if(0==nasline.find('*')):
            #iPt = int(nasline[1:8])
            mapResults['pts'][iPt-1,2]=float(nasline[8:24])
        else:
            for type in mapNAS2NP.keys():
                elems_strs = mapNAS2NP[type][1]
                vertex_num = mapNAS2NP[type][0]
                if(0==nasline.find(type)):
                    iElem = mapResults['i%s'%elems_strs]
                    # print('iElem%d'%iElem)
                    mapResults['i%s'%elems_strs] = iElem + 1
                    mapResults[elems_strs] = double_size(mapResults[elems_strs],iElem+1,1.0)
                    mapResults[elems_strs][iElem,vertex_num] = int(nasline[16:24])
                    # print(vertex_num)
                    for iVertex in np.arange(0,vertex_num):
                        beg = 8*(int(iVertex)+3)  #24
                        end = 8*(int(iVertex)+4)  #24
                        # print(iVertex,beg,end)
                        # print('%s %d:%d'%(type,beg,end))
                        mapResults[elems_strs][iElem,iVertex] = int(nasline[beg:end])-1
        if(0==nasline.find('ENDDATA')):
            break
    mapResults['pts'] = mapResults['pts'][:max_ipt,:]
    for type in mapNAS2NP.keys():
        elems_strs = mapNAS2NP[type][1]
        vertex_num = mapNAS2NP[type][0]
        nElems = mapResults['i%s'%elems_strs]
        # print (type,nElems)
        mapResults[elems_strs] = mapResults[elems_strs][:nElems,:]
    return mapResults
def fileExt(file_path):
    import os
    # 使用 splitext 分割文件名和后缀
    _, ext = os.path.splitext(file_path)
    # 返回后缀（不带点，例如 "txt"）
    return ext.lstrip('.') if ext else ''
def console_cmd(cmd,show_on_console=True,wait=False):
    import subprocess
    file = 'my_console_out.log'
    if(show_on_console):
        file=None
    else:
        file=open(file,'a')
        file.write(cmd+"\n")
    status = subprocess.Popen(cmd,bufsize=0,stdout=file)
    if(show_on_console==False):
        file.close()
    if(wait==True):
        status.wait()
def line2strs(line,prefix=''):
    line = line.replace(prefix,'')
    #每一行字符串的前后各加一个空格
    line=' '+line+' '
    #每一行字符串的'(',')','      '用' '替换
    sepList = ['\n','\t','(',')',',','     ','    ','   ','  ',]
    for sep in sepList:
        line = line.replace(sep,' ')
    strs = line.split(' ')[1:-1]
    return strs
def line2data(line,prefix=''):
    '''
    将一行数据分割成字符串列表
    '''
    strs=line2strs(line,prefix)
    return np.array(list(map(lambda u:float(u),strs)))
def line1stData(line):
    '''
    读取某行字符串中指定分割符的第一个数值
    '''
    strs = line2strs(line)
    return float(strs[0])
def lines2data(lines,cnt=-1):
    if cnt==-1:
        cnt = int(line1stData(lines[0]))
    results = np.zeros((0,3))
    for ind in np.arange(cnt):
        if ind==0:
            column = np.size(line2data(lines[ind+1]))
            results = np.zeros((cnt,column))
        results[ind,:]=line2data(lines[ind+1])
    return results
def read_facets(path):
    '''
    从wires文件中提取points和wires
    '''
    lines = file2lines(path)
    pt_cnt = int(line1stData( lines[4]))
    pts = lines2data(lines[4:(5+pt_cnt)])
    lines = lines[7+pt_cnt:]
    strs = line2strs(lines[0])
    nprop=(int(strs[1]))
    wire_cnt = int(strs[0])
    wires = lines2data(lines)[:,:nprop+1]
    wires[:,nprop:]=wires[:,nprop:]+1
    return (pts,wires.astype(np.int64)-1)

def mesh_filter_sphere(mesh_file,rx,radius):
    '''
    以rx为中心，以radius为半径的球面来截取mesh_file中的网格，进行提取目标点
    返回球内的目标点
    '''
    ext = fileExt(mesh_file)
    if(ext=="nas"):
        mapResults = read_nas(mesh_file)
        tris=mapResults['tris']
        pts = mapResults['pts']
    elif(ext=="facet"):
        pts,tris = read_facets(mesh_file)
    x0,y0,z0 = rx
    ipts_filter = set()
    for tri in tris:
        tri_ctr = 1.0/3*(pts[tri[0],:3] + pts[tri[1],:3] + pts[tri[2],:3])
        r0 = np.sqrt((tri_ctr[0]-x0)**2.0+(tri_ctr[1]-y0)**2.0+(tri_ctr[2]-z0)**2.0)
        if(r0<radius):
            ipts_filter.add(tri[0])
            ipts_filter.add(tri[1])
            ipts_filter.add(tri[2])
    npts_filter = len(ipts_filter)
    pts_filter = np.zeros((npts_filter,3))
    # for ipts in ipts_filter:
    #     pts_filter[]
    pts_filter = pts[list(ipts_filter),:]
    return pts_filter

def tri6_N(ξ, η):
    return np.array([ξ * (2 * ξ - 1),η * (2 * η - 1),2 * (1 - ξ - η)**2.0 - (1 - ξ - η),4 * ξ * η,4 * η * (1 - ξ - η),4 * ξ * (1 - ξ - η)])

def tri6_interp(vertex, ξ, η):
    bas = tri6_N(ξ, η)
    result = np.array([0.0,0.0,0.0])
    for ibas in np.arange(6):
        result += bas[ibas] * vertex[ibas,:]
    return result